import math


class Converter:
    def __init__(self) -> None:
        self.coord

    def coord(self, coord):
        """
        实例对象中传入坐标
        @params coord：[x,y,z] 三坐标值
        @return self
        """
        self.coord = coord
        return self

    def geodetic_spatialRectangular(self, crs):
        """
        大地坐标系转空间直角坐标系
        @parmas crs：crs() 椭球参数，详见crs方法
        @return self
        """
        longitude, latitude, height = (
            math.radians(float(self.coord[0])),
            math.radians(float(self.coord[1])),
            float(self.coord[2]),
        )

        E1S = crs["FIRST_ECCENTRICITY_SQUARE"]  # 椭球第一偏心率的平方
        a = crs["SEMIMAJOR_AXIS"]  # 长半轴
        W = math.sqrt(1 - E1S * math.sin(latitude) * math.sin(latitude))
        N = a / W  # 卯酉圈半径  N=a/(1-e²sin²B)^(1/2)

        # 计算 X Y Z
        X = (N + height) * math.cos(latitude) * math.cos(longitude)  # X=(N+H)cosBcosL
        Y = (N + height) * math.cos(latitude) * math.sin(longitude)  # Y=(N+H)cosBsinL
        Z = (N * (1 - E1S) + height) * math.sin(
            latitude
        )  # Z=(N(1-e²)+H)sinB  -----e 椭球第一偏心率

        self.coord = [X, Y, Z]
        return self

    def spatialRectangular_geodetic(self, crs):
        """
        空间直角坐标系转大地坐标系
        @parmas crs：crs() 椭球参数，详见crs方法
        @return self
        """
        X, Y, Z = self.coord[0], self.coord[1], self.coord[2]
        E1S = crs["FIRST_ECCENTRICITY_SQUARE"]  # 椭球第一偏心率的平方
        a = crs["SEMIMAJOR_AXIS"]  # 椭球长半轴

        # 初始纬度（弧度）
        tan_latitude_init = Z / math.sqrt(
            X * X + Y * Y
        )  # 迭代初始纬度 tanB = Z / (X² + Y²)^(1/2)
        # 迭代计算
        while True:
            N_init = a / math.sqrt(
                1
                - E1S
                * math.sin(math.atan(tan_latitude_init))
                * math.sin(math.atan(tan_latitude_init))
            )
            tan_latitude_final = (
                Z + N_init * E1S * math.sin(math.atan(tan_latitude_init))
            ) / (math.sqrt(X * X + Y * Y))
            if abs(tan_latitude_final - tan_latitude_init) <= 1e-20:
                break
            tan_latitude_init = tan_latitude_final
        latitude = math.atan(tan_latitude_init)

        # 经度（弧度）
        longitude = math.atan(Y / X)
        # 高度（米）
        W = math.sqrt(1 - E1S * math.pow(math.sin(latitude), 2))
        N = a / W
        height = math.sqrt(X * X + Y * Y) / math.cos(latitude) - N
        # 弧度转角度
        longitude = math.degrees(longitude)
        if longitude < 0:
            longitude = longitude + 180
        latitude = math.degrees(latitude)

        self.coord = [longitude, latitude, height]
        return self

    def bursa7params_converte(self, bursa_params):
        """
        不同椭球的空间直角坐标系 布尔莎其参数转换
        @params bursa_params：bursa_params()
        return self
        """
        x, y, z = self.coord[0], self.coord[1], self.coord[2]
        dx, dy, dz, rx, ry, rz, k = (
            bursa_params[0],
            bursa_params[1],
            bursa_params[2],
            math.radians(bursa_params[3] / 3600),
            math.radians(bursa_params[4] / 3600),
            math.radians(bursa_params[5] / 3600),
            bursa_params[6] * (1e-6),
        )

        a1 = k + 1
        a2 = a1 * rx
        a3 = a1 * ry
        a4 = a1 * rz

        x_r = dx + a1 * x - a3 * z + a4 * y
        y_r = dy + a1 * y + a2 * z - a4 * x
        z_r = dz + a1 * z - a2 * y + a3 * x
        self.coord = [x_r, y_r, z_r]
        return self

    def positive_gaussain_projecton(self, crs):
        """
        高斯正算(3度带投影)
        @parmas crs：crs() 椭球参数，详见crs方法
        @return self
        """
        a = crs["SEMIMAJOR_AXIS"]
        E1S = crs["FIRST_ECCENTRICITY_SQUARE"]

        longitude, latitude, height = self.coord[0], self.coord[1], self.coord[2]

        latitude_rad = math.radians(latitude)

        zone_num = int((longitude + 1.5) / 3.0)  # 计算3度带投影度带号
        central_meridian = zone_num * 3  # 计算中央经线
        dLongitude = longitude - central_meridian  # 经差
        tsin = math.sin(latitude_rad)
        tcos = math.cos(latitude_rad)
        t = math.tan(latitude_rad)
        m = (math.pi / 180.0) * dLongitude * tcos
        et2 = E1S * pow(tcos, 2)
        et3 = E1S * pow(tsin, 2)
        X = (
            111132.9558 * latitude
            - 16038.6496 * math.sin(2 * latitude_rad)
            + 16.8607 * math.sin(4 * latitude_rad)
            - 0.0220 * math.sin(6 * latitude_rad)
        )
        N = a / math.sqrt(1 - et3)

        x = X + N * t * (
            0.5 * pow(m, 2)
            + (5.0 - pow(t, 2) + 9.0 * et2 + 4 * pow(et2, 2)) * pow(m, 4) / 24.0
            + (61.0 - 58.0 * pow(t, 2) + pow(t, 4)) * pow(m, 6) / 720.0
        )
        y = 500000 + N * (
            m
            + (1.0 - pow(t, 2) + et2) * pow(m, 3) / 6.0
            + (5.0 - 18.0 * pow(t, 2) + pow(t, 4) + 14.0 * et2 - 58.0 * et2 * pow(t, 2))
            * pow(m, 5)
            / 120.0
        )
        self.coord = [x, y, height]
        return self

    def negative_gaussian_projection(self, crs, central_meridian=120):
        """
        高斯反算(3度带投影)
        @paramas crs：crs() 椭球参数，详见crs方法
        @params central_meridian：int 中央子午线 默认120°E
        @return self
        """
        X = self.coord[0]
        Y = self.coord[1]
        height = self.coord[2]
        a = crs["SEMIMAJOR_AXIS"]
        E1S = crs["FIRST_ECCENTRICITY_SQUARE"]  # 第一偏心率平方
        EE = (1.0 - math.sqrt(1 - E1S)) / (1.0 + math.sqrt(1 - E1S))
        E2S = E1S / (1 - E1S)  # 第二偏心率平方

        ProjNo = int(X / 1000000)
        L0 = math.radians(central_meridian)
        X0 = ProjNo * 1000000 + 500000
        Y0 = 0
        xval = X - X0
        yval = Y - Y0

        M = yval
        u = M / (a * (1 - E1S / 4 - 3 * E1S * E1S / 64 - 5 * E1S * E1S * E1S / 256))

        # 偏导数
        fai = (
            u
            + (3 * EE / 2 - 27 * EE * EE * EE / 32) * math.sin(2 * u)
            + (21 * EE * EE / 16 - 55 * EE * EE * EE * EE / 32) * math.sin(4 * u)
            + (151 * EE * EE * EE / 96) * math.sin(6 * u)
            + (1097 * EE * EE * EE * EE / 512) * math.sin(8 * u)
        )
        C = E2S * math.cos(fai) * math.cos(fai)
        T = math.tan(fai) * math.tan(fai)
        NN = a / math.sqrt(1.0 - E1S * math.sin(fai) * math.sin(fai))
        R = (
            a
            * (1 - E1S)
            / math.sqrt(
                (1 - E1S * math.sin(fai) * math.sin(fai))
                * (1 - E1S * math.sin(fai) * math.sin(fai))
                * (1 - E1S * math.sin(fai) * math.sin(fai))
            )
        )
        D = xval / NN

        # 计算经纬度（弧度）
        longitude_rad = L0 + (
            D
            - (1 + 2 * T + C) * D * D * D / 6
            + (5 - 2 * C + 28 * T - 3 * C * C + 8 * E2S + 24 * T * T)
            * D
            * D
            * D
            * D
            * D
            / 120
        ) / math.cos(fai)
        latitude_rad = fai - (NN * math.tan(fai) / R) * (
            D * D / 2
            - (5 + 3 * T + 10 * C - 4 * C * C - 9 * E2S) * D * D * D * D / 24
            + (61 + 90 * T + 298 * C + 45 * T * T - 256 * E2S - 3 * C * C)
            * D
            * D
            * D
            * D
            * D
            * D
            / 720
        )

        # 弧度转角度
        longitude = math.degrees(longitude_rad)
        latitude = math.degrees(latitude_rad)

        self.coord = [longitude, latitude, height]
        return self

    def export_Coordinate(self, length=[8, 8, 4]):
        """
        输出类属性中存储的坐标
        @params length: int[len,len,len] 三坐标小数点保留的长度
        @return self.coord： int[x,y,z] 三坐标
        """
        return [
            round(self.coord[0], length[0]),
            round(self.coord[1], length[1]),
            round(self.coord[2], length[2]),
        ]

class CRS_Builder:
    WGS84 = {
        "SEMIMAJOR_AXIS": 6378137.0,
        "OBLATENESS_RECIPROCAL": 298.257223563,
        "FIRST_ECCENTRICITY_SQUARE": 0.00669437999013,
    }

    CGCS2000 = {
        "SEMIMAJOR_AXIS": 6378137.0,
        "OBLATENESS_RECIPROCAL": 298.257222101,
        "FIRST_ECCENTRICITY_SQUARE": 0.00669438002290,
    }

    XIAN80 = {
        "SEMIMAJOR_AXIS": 6378140.0,
        "OBLATENESS_RECIPROCAL": 298.257,
        "FIRST_ECCENTRICITY_SQUARE": 0.00669438499959,
    }


class Bursa_Params_Builder:
    WGS84_CGCS2000 = [
        178.868103027344,
        96.256584167805,
        5.05481958389282,
        -1.684915,
        3.413315,
        -5.365031,
        -0.386750966310501,
    ]
    WGS84_XIAN80 = [
        152.225021,
        94.981697,
        5.499349,
        -0.748785,
        1.118085,
        -1.486786,
        -1.486786,
    ]
